Mon. Not. R. Astron. Soc. OOO.HlfUl (20101 Printed 3 March 2013 (MN KTriX style file v2.2) 



ROPS: A New Search for Habitable Earths in the Southern Sky 

J.R. Barnes 1 , J.S. Jenkins 2 , H.R.A. Jones 1 , P. Rojo 2 , R Arriagada 3 , A. Jordan 3 , 
- - - . D. Minniti 3 ' 4 , M. Tuomi 15 , S.V. Jeffers 6 , and D. Pinfield 1 

N i 1 Centre for Astrophysics Research,. University of Hertfordshire,. College Lane, Hatfield. Herts. ALIO 9AB. UK. 
' - _J ' 2 Departamento de Astronomi'a, Universidad de Chile, Camino del Observatorio 1515, Las Condes, Santiago. Chile. 
^■"^ i 3 Dept. of Astronomy and Astrophysics, Pontificia Universidad Catolica de Chile, Av. Vicuna Mackenna 4860, Santiago, Chile 
V N i 4 y a fi can Observatory, V00120 Vatican City State, Italy 

5 University of Turku, Tuorla Observatory, Department of Physics and Astronomy, Vdisdldntie 20, FI-21500, Piikkid, Finland 
' 6 Institut fur Astrophysik, Georg-August-Universitdt, Friedrich-Hund-Platz 1, Friedrich-Hund-Platz 1, D-37077 Gttingen. Germany. 



Submitted March 2012; Accepted May 2012 



Oh 



ABSTRACT 

We present the first results from our Red Optical Planet Survey (ROPS) to search for low mass 
' planets orbiting late type dwarfs (M5.5V - M9V) in their habitable zones (HZ). Our observa- 

tions, with the red arm of the MIKE spectrograph (0.5 - 0.9 |J.m) at the 6.5 m Magellan Clay 
telescope at Las Campanas Observatory indicate that ^ 92% of the flux lies beyond 0.7 \xm. 
£^ ' We use a novel approach that is essentially a hybrid of the simultaneous iodine and ThAr 

£^ , methods for determining precision radial velocities. We apply least squares deconvolution to 

■ obtain a single high S/N ratio stellar line for each spectrum and cross correlate against the 

simultaneously observed telluric line profile, which we derive in the same way. 

Utilising the 0.62 - 0.90 [im region, we have achieved an r.m.s. precision of 10 ms _1 for 
an M5.5V spectral type star with spectral S/N ~ 160 on 5 minute timescales. By M8V spec- 
tral type, a precision of ~ 30 ms _1 at S/N = 25 is suggested, although more observations 
qq ■ are needed. An assessment of our errors and scatter in the radial velocity points hints at the 

' presence of stellar radial velocity variations. Of our sample of 7 stars, 2 show radial veloc- 

\Q , ity signals at 6a and 10er of the cross correlation uncertainties. We find that chromospheric 

activity (via Ha variation) does not have an impact on our measurements and are unable to 
f~>) | determine a relationship between the derived photospheric line profile morphology and radial 

. velocity variations without further observations. If the signals are planetary in origin, our find- 

ings are consistent with estimates of Neptune mass planets that predict a frequency of 13 - 27 
per cent for early M dwarfs. 

Our current analysis indicates the we can achieve a sensitivity that is equivalent to the 
K> ' amplitude induced by a 6 M® planet orbiting in the habitable zone. Based on simulations, we 

5_j , estimate that <10 M® habitable zone planets will be detected in a new stellar mass regime, 

with ^20 epochs of observations. Higher resolution and greater instrument stability indicate 
that photon limited precisions of 2 ms" 1 are attainable on moderately rotating M dwarfs (with 
v sin i ^5 kms -1 ) using our technique. 

Key words: (stars:) planetary systems stars: activity stars: atmospheres stars: spots tech- 
niques: radial velocities 
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1 INTRODUCTION 

It is reasonable to expect that planets with predominantly lower 
mass may be found orbiting low mass stars. Core accretion theory 
predicts that the formation o f Jupiter mass planet s orbiting M dwarf 
stars is seriously inhibited dLaughlin et al.ll2004b . This argument is 
in good agreement with the Keck survey observations which esti- 
mate that the fraction of Jupiter-mass plane ts orbiting M, K an d 
G stars is 1.8%, 4.2% and 8.9% respectively djohnson et al.ll2007b . 
More recently, it has been shown that the occurrence rate of gas 
giants orbits of less than 2000 days is 3-10 times smaller for M 



dwarfs than for F, G and K dwarfs dCumming et~ail 2008). In fact, 
models predict even lower incidences of Jupiter-mass planets, al- 
though the trend of a decrease in numbers with decreasing s tellar 
mass agrees with observations. For instance. Ilda & Liril d2005h pre- 
dict ~ 7.5 times lower incidence of hot Jupiters around 0.4 M© 
when compared with 1.0 M© and lKennedv & Kenvod 120081 ) pre- 
dict a factor of ~ 6 difference. 

However, fo r lower mass protoplanetary disks, predictions by 
Ida & Liril J20051) indicate that although ice cores can form, they can 
still acquire mass through inefficient gas accretion. The resulting 
bodies form Neptune mass (or ~ 10 Mq) planets which are able 
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to undergo type II migration, ending in close proximity to the par- 
ent star, with radii < 0.05 AU. It has also been noted that the semi- 
major axis distribution of known radial velocity planets o rbiting 
M sta rs peaks at closer orbits than for more massive stars dCurriel 
2009). While it is expected that the frequency of Neptune mass 
close-in planets peaks for M ~ 0.4 Mq stars, the mass distribution 
peak is expected at 10 M© for 0.2 Mq stars, with a tail extending 
down to a few M®. In fact the Kepler planetary candidates now 
confirm that Meto Mjv ep t planets are ~ 2.5 - 4 times more numer - 
ous around M stars than earlier spectral types dBorucki et alj201lh . 
and crucially are found in greater numbers than transiting Jupiter 
size planets orbiting earlier spectral types. Many of these candi- 
dates are in very close orbits of < 0. 1 AU. Moreover, there are in- 
dications faoward et alj2010Mwittenmver et al.l201 ll) that 1 1.8%- 
17.4% of solar-like stars possess rocky (1-10 Msarth) planets. 
Scaling these values by the 2.5 - 4 times increased rocky planet oc- 
currence rate reported by Kepler indicates that 30 - 7 per cent of 
M dwa rfs should be orbited by close-in rocky planets. Bon fils et al.l 
d201 ll) have recently reported on occurrence rates from the early- 
M dwarf HARPS (High Accuracy Radial Velocity Planet Searcher) 
sample and find a similarly high occurrence rate for rocky planets 
(77©) of 0.5 4 j^. Similarly, for Neptune mass planets, s caling the 
respective (How ard et al. I I2OIOI: IWittenmver et al.l l201lh frequen- 
cies of 6.5 and 8.9 per cent by the Kepler candidate planet frequen- 
cies leads to expected occurrence rates of 13-27 per cent for M 
dwarfs. Whether rocky planets that orbit in the classical habitable 
zone (where liquid water could be present) are in reality habitable, 
may depend on fac t ors such as tidal locking and stellar activity 
jKastingetaT]ll99l l Tarte r et al. More recent studies, that 

investigate factors such as obliq uity may al so be important consid- 
erations for habitability dHeller et al.ll20l"lh . although these topics 
are beyond the scope of this paper which focusses merely on plan- 
etary detection. 

Despite comprising 70% of the solar neighbourhood popula- 
tion, the lower mass M dwarfs have remained beyond efficient de- 
tection with optical based surveys as the stellar flux peaks at longer 
wavelengths. At infrared wavelengths, surveys targeted at search- 
ing for low mass planetary systems associated with low mass stars 
can be achieved with realistic ti mescales. Neverth eless, to date, 
only one survey, using CRIRES teeanetalj|2010l) . has reported 
sub- 10 ms" 1 precision, with 5.4 ms" 1 reported using an NH3 gas 
cell and modest 364 A of spectrum in the K-band. More recent 
surveys have used telluric lines as a reference fiducial, and while 
iRodler et alj fe012h hav e obta ined 180-300 ms 1 precision at R 
-20,0000, Bailey etal.|j2012h have achieved ~50 ms 1 precision 
which is limited by activity on their sample of young active stars. 
Equally, the red-optical regime offers similar advantages with a 
V - 1 ~ 2 - 5 for the earliest-latest M dwarfs respectively, but with 
m uch improved wav elength coverage over the near-infrared regime 
of lBeanetai1 ( l2010h . There are also significant telluric free spectral 
windows in the red-optical. 

The effects of jitter due to stellar activity are expected to 
present limitations but are also expected to be lower for M 
dwarfs than for F, G & K dwarfs since the contrast ratio between 
spots and stellar photosphere is lower and the starspot distribu- 
tions may be more uniformly dist ributed across the stellar surface 
jBarnes & Collier Camero nl boOll) . We have simulated the radial 
velocity (RV) noise for M dwarfs (using a 3D, 2-temperature stellar 
model) with low (solar) activity levels and estimate that the radial 
velocity jitter is ~ 1 .5 - 2 times smaller in the I-band vs the V-band 
jBarnes et al.ll201 il) . We estimate that the appropriate jitter in the 
near infrared Y-band for even an extreme solar activity M dwarf 
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Figure 1. The mean flux level of the 33 extracted MIKE orders cover- 
ing 4830-9172 A for M5.5V and M9V targets (scaled to the same level 
at I band centre). Also shown are the optical I2 cell regime and the 
HARPS coverage that overlap with the MIKE red arm data. We find that 
F0.7-O.9/F0.5— 0.7 = 1 1.5 and 19 for the M5.5V and M9V spectra respec- 
tively, demonstrating the inefficiency of wavelengths < 7000 Afor observ- 
ing M dwarf targets. 



analogue with v sin i = 5 kms is in the range ~0.7 - 5 ms -1 for 
an M6V star. The uncertainty arises from estimates of extremes of 
spot contrast, with 0.7 ms -1 corresponding to a photosphere-spot 
temperature difference of T p hot - T spot = 200 K, while 5 ms 
arises when the contrast is much higher, with T spo t/T p hot = 0.65. 
Given th at mean solar spot tem peratures are of order T p hot - T spo t 
= 500 K (Lagrange et alj20ld ). we can assume that ~2 - 3 ms" 1 is 
a reasonable estimate of the spot induced jitter of a 5 kms" 1 star 
when we scale our results to the I-band. 

Flaring is known to be an important contributor to variations 
in spectral features in stars, but being high energy phenomena tend 
to affect bluer wavelen gths more than redder and infrared wave- 
lengths. iReinerj d2009h investigated flaring in the nearby active 
M6 dwarf, CN Leo, from nearly 200 observations spanning three 
nights. The spectra, taken with uves at the Very Large Telescope 
array were taken at at R = 40,000, covering 0.64- 1.08 u.m, and 
thus closely resemble our observation setup with MIKE (f|2j. It 
was found that a large flaring event, easily identifiable by observ- 
ing strong transient emission in Ha, could effect radial velocity 
measurements by ~500 ms" 1 in orders containing strong emis- 
sion features. After the flaring event, the precision was limited to 
~50 ms" 1 . In other orders however, it was found that the radial ve- 
locity jitter is below the 10 ms" 1 level, even during flaring events. 
Careful choice of line regions is therefore important and inspection 
of activity indicators, such as Ha important where sub-10 ms -1 
precision is required. 

In section ij2]we discuss our observations and data reduction 
and demonstrate the case for attempting precision radial velocities 
of M dwarfs at red-optical wavelengths. In ij3] we justify our choice 
of reference fiducial and give details of our wavelength calibration 
procedure. Sj4] introduces the least squares deconvolution approach 
to measuring precision radial velocities before we present results 
from simulations and our observations in Sj5] In light of our obser- 
vational radial velocities, a discussion (ij6]l that further considers 
the sources of noise in our results is given, before we make make 
concluding statements of the prospects of this kind of survey (f|7J. 
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2 OBSERVATIONS 

We observed a number of bright M dwarf targets with the Mag- 
ellan Inamori Kyocera Echelle (MIKE) spectrograph at the 6.5 m 
Magellan Clay telescope on 2010 November 21 & 22. Although 
the instrument simultaneously records two wavelength regions in a 
blue arm 3350 - 5000 A and a red arm 4900 - 9500 A, we only made 
use of the latter owing to the low S/N ratio obtained on our mid-late 
M dwarf targets. We observed with a 0.7" slit. 

Since ThAr lamps exhibit many lines for calibration, and are 
generally always used with echelle spectrographs working at opti- 
cal wavelengths, we made regular observations with the compari- 
son lamp avai l able w ith MIKE. HARPS spectra have been used by 
lLovis & Pepel d20070 to determine ThAr line wavelengths with a 
precision of ~ 10 ms ^(O.IS mA) on average for individual lines. 
This enables global wavelength solutions that attain sub-ms -1 sta- 
bility when several hundred lines, spanning many echelle orders, 
are used. The median FWHM of the ThAr lines was found to be 
4.05 pixels, indicating a mean resolution of R = 34,850. The ob- 
serving conditions over the two nights were very good, with seeing 
estimates in the range 0.7 - 1.2 for targets observed at airmasses < 
1.5. Our targets are listed in Table [T] and with the exception of 
LP944-20 are all known to exhibit vsini ^ 5 kms -1 (Jenkins et 
al., In prep). Our observing sequence comprised of taking spectra 
alternated with calibration frames since it is not possible obtain a 
simultaneous reference observation with MIKE at red wavelengths. 
We note that MIKE does possess an iodine cell but that the wave- 
length cutoff of ~ 7000 A means that we are unable to make use of 
I2 lines. 

2.1 Data extraction 

Pixel to pixel variations were corrected for each frame using flat- 
field exposures taken with an internal tungsten reference lamp. 
Since few counts are recorded in the reddest orders of the red CCD 
on MIKE when care is taken to keep counts in the middle orders 
to no more than 50,000 counts (full well capacity 65,536 counts), 
we resorted to taking two sets of flat field frames of different ex- 
posure lengths. A total of 60 exposures of 40s each were used to 
flatfield the reddest orders. The worst cosmic ray events were re- 
moved at the pre-extraction sta ge using the Starlink FIGARO rou- 
tine BCLEAN dShortridgell 1 9931) . The spectra were extracted using 
ECHOMOP's implementatio n of the optimal extraction algorithm 
developed by Home ( 1986). ECHOMOP rejects all but the strongest 
sky lines I IBarnes et al.ll2007bl) and propagates error information 
based on photon statistics and readout noise throughout the extrac- 
tion process. 

2.2 M dwarf spectral fluxes at red-optical wavelengths 

Fig.Qjcompares the wavelength regimes used to make precision ra- 
dial velocities. We note that by utilising the 0.62 - 0.90 p.m region, 
we use almost the same wavelength extent as HARPS which covers 
0.38-0.69 pin. The fluxes observed with MIKE are shown for an 
M5.5V spectrum and an M9V spectrum. The spectral energy dis- 
tributions are determined from the observations with no flux cali- 
bration, and thus include the spectral energy distribution of the stars 
and the response of the telescope and instrument (i.e. CCD quan- 
tum efficiency and instrumental efficiency and blaze function). The 
curves are derived by taking the mean counts in each spectral or- 
der recorded on the red arm CCD of MIKE. The benefit of using 
red-optical wavelengths for mid-late M dwarfs is obvious with a 



recorded flux in the 0.7 - 0.9 pm regime that exceeds the flux mea- 
sured at 0.5 - 0.7 (xm (traditionally used by surveys that utilise ThAr 
or iodine reference lines), by a factor of 1 1.5 (M5.5V) to 19 (M9V). 
This is equivalent to a magnitude advantage of ~2.65 - 3.2. 

The Doppler information available in these regimes is equally 
important and must be taken into consideration when estimating 
velocity precisions achievable in different wavelength regimes, as 
found by iReiners et al. for example. For a M5.5V dwarf, 

the total number of lines available to HARPS is ~ 16,000, while the 
0.6-0.9 pm region contains ~ 12,300 lines that are on average 14 
per cent less deep dBrott & Haus childt 2005). Similarly, the num- 
ber of lines in the 0.5 - 0.7 [im region plotted in Fig.Qjis ~ 10,3000, 
while only ~ 7400 lines are available in the 0.6 - 0.9 ixm region, 
with relative normalised depths of 0.67. For an M5.5V dwarf, as- 
suming Poisson statistics correlate linearly with precision, the rel- 
ative precision for a fixed integration time in the red-optical is ex- 
pected to be 2.4 times that attained at standard optical wavelengths 
with the same integration time. In other words, to obtain the same 
precision at standard optical wavelengths would require 5.6 times 
the integration time . Simi larly, for an M9V dwarf, the models of 
iBrott & Hauschildtl {2005) indicate line number and depth ratios 
(0.7-0.9 ixm / 0.5-0.7 urn) of 0.755 and 0.925 respectively. Pre- 
cision increases by 3.6 times when working in the red-optical, or 
equivalently integration times of only ~ 1/13 are required to achieve 
the same precision as standard optical wavelength surveys. 



3 WAVELENGTH CALIBRATION 

Wavelength calibrations were made by an initial order-by-order 
identificat ion of suitable lines using the ThAr wavelengths pub- 
lished by lLovis & Pepel J2007h . which are estimated to enable a 
calibration (i.e. global) r.m.s to better than 20 cms -1 for HARPS 
(see §2). Pixel positions were initially identified for a single arc us- 
ing a simple Gaussian fit. For each subsequent arc, a cross-match 
was made followed by a multiple-Gaussian (up to three profiles) fit 
around e ach identified line using a Levenberg-Marquardt fitting al- 
gorithm dPress et alJI 19860 to obtain the pixel position of each line 
centre. The Lovis & Pepe line list was optimised for HARPS at R 
= 110,000, while our observations were made at R ~ 35,000 ne- 
cessitating rejection of some lines. Using a multiple Gaussian fit 
enables the effect of any nearby lines to be accounted for in the 
fit that also included a first order (straight line) background. Any 
lines closer than the instrumental FWHM were not used. Finally 
for each order, any remaining outliers were removed after fitting a 
cubic-polynomial. In addition, any lines that were not consistently 
yielding a good fit for all arc frames throughout both nights (to 
within 3-sigma of the cubic fits) were removed. A total of 32-71 
lines were available for the orders centered at 6363 - 7470 A. For 
the redder orders centered at 7636 - 88 1 1 A, 1 8 - 25 lines were avail- 
able. To improve stability, we made a two dimensional wavelength 
calibration, with 4 coefficients in the dispersion direction and 7 co- 
efficients in the cross-dispersion direction. By iteratively rejecting 
outlying pixels from the fit, we found that of the input 575 lines, 
clipping the 15 furthest outliers yielded the most consistent fit from 
one solution to the next. Sigma clipping has the potential to remove 
blocks of lines in a given region and we found this to give less sta- 
ble solutions. Of the 575 input lines, we clipped 15 from each frame 
leaving a total of 550 lines for each solution. The zero point r.m.s. 
(i.e. the r.m.s. by combining all lines) is found to be 8.73 ± 0.32 
ms - . Wavelength calibration precision could be improved, in the 
reddest orders at least, where A > 0.85 pm, by simultaneously us- 
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Star 


SpT 


Imag 


Exp 
[s] 


v sin i 
kms 


S/N (extracted) 


S/N (deconvolved) 


Number of observations 


GJ 1002 


M5.5V 


10.2 


300 


^5 


113 


4610 


8 


GJ 1061 


M5.5V 


9.5 


300 


<5 


163 


7400 


12 


GJ 1286 


M5.5V 


11.1 


350 


<5 


84 


3480 


8 


GJ3128 


M6V 


11.1 


350 


<5 


84 


3480 


8 


SO J025300.5+165258 


M7V 


10.7 


350 


<5 


90 


3150 


9 


LHS 132 


M8V 


13.8 


500 


<5 


25 


875 


4 


LP944-20 


M9V 


13.3 


500 


30 


44 


1540 


7 



Table 1. List of targets observed on November 21 & 22 with their estimated spectral types, I band magnitudes, v sini values (Jenkins et al., In prep.) and 
exposure times. Extracted S/N ratio and S/N ratio after deconvolution are tabulated in columns 6 & 7. Column 8 lists the total number of observations on each 
target. 



ing a ThAr and UNe lamps. The UNe lamp has been shown to pro- 
vide around 6 times more lines in this region than the ThAr lamp 
dRedmanetal]201lh . 

We find that despite a precision of sub 10 ms" 1 on a given 
wavelength solution, the solution from one ThAr frame to the next 
yields r.m.s. residuals of 30-80 ms -1 . The difference does not 
show a constant shift or tilt across all orders but may appear random 
from one order to the next indicating that the long term solution is 
in fact not stable. The dynamic range of the ThAr lines used was 
however large, so that weak lines, especially in the redder order 
where fewer lines per order are available, contribute to a less stable 
solution that is desirable. A total of 88 per cent of the lines possess 
strengths that are ^ 0.05 of the maximum line used (even consider- 
ing the removal of lines that peaked with more than 50000 counts 
before extraction), with most of these lines possessing only a few 
thousand counts above the bias level. 



3.1 Stability of mike 

In Fig. [2] we illustrate the stability of MIKE by cross-correlating 
ThAr frames. This enables us to ascertain our precision when in- 
terpolating our reference velocity for each target observation. The 
MIKE User's Guid^] indicates that the instrument is stable at the 
~ 1 pixel level (in the cross-dispersion direction), owing to small 
temperature changes through the night. However, Fig. [2]in fact in- 
dicates that use of MIKE as a precision radial velocity instrument 
requires that a simultaneous reference fiducial be used if consis- 
tent results are to be achieved. The crosses joined by solid lines 
represent the drift in ms -1 as a function of frame number (cov- 
ering all science frames taken throughout the night). The hatched 
regions represent times at which the telescope was pointing at the 
same right ascension and declination; in other words, no slewing 
occurred during these intervals. Nevertheless, it can be seen that 
large shifts of the instrument take place during these intervals. In 
some cases, particularly on the second night, the wavelength ap- 
pears to have drifted by ~ 500 ms -1 , or ~ 0.25 pixels between arc 
observations. In many cases, the drift is somewhat less, but clearly 
severely prohibits our ability to make precision radial velocity mea- 
surements of ^ 10 ms -1 by relying on ThAr measurements. 

Moreover, inspection of telluric lines indicates changes in po- 
sition that do no correlate with the arc lines. We suspect that this 
apparent discrepancy, and the semi-random nature of the measured 
offsets may arise from mechanical movement of the calibration 
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mirror when alternating between ThAr observations and science 
target observations. 



3.2 Precision velocity measurements with mike 

While Fig. [^indicates that the ThAr lines can shift by large quan- 
tities, Fig.[3]demonstrates the relationship between the ThAr shifts 
and the shifts from telluric lines (i.e. the telluric lines do not show 
the same shifts as the ThAr lines). The plot shows the relative drift 
in pixels as a function of pixel for all orders, clearly indicating that 
a tilt is also present. We measured the strongest telluric lines from 
observations of GJ1061 since this object has the highest S/N ~ 160. 
We measured the telluric line positions in the same manner as de- 
scribed above for the ThAr frames. Plotted in Fig. [3]are the shifts 
relative to the first observation of GJ 1061 on the first night. It can 
be seen that there is a significant shift, even between the pair of 
observations, which were made together with a ThAr observation 
between them. Similarly, we compared the ThAr frame between 
these two observations with that taken immediately after the first 
GJ 1061 observation. The scatter in the telluric lines is too great to 
give a reliable tilt measurement, but the ThAr relationship is much 
tighter. Nevertheless, the tilt appears to be approximately the same 
for the tellurics and the ThAr frames, even of the shift is different. 
The similarity of tilt, but differing shift is seen for other pairs of ob- 
servations of GJ 1061. However, owing to the larger scatter seen in 
telluric pixel positions for all other objects (which typically possess 
0.5 - 0.75 of the S/N of our GJ1061 observations), we are unable to 
easily see the tilt. 

Our strategy for making precision radial velocities stems from 
the above observations. We assume that the local wavelength solu- 
tion can be described as a shift and a tilt relative to a single master 
spectrum. Since we find that the local 2D wavelength solution re- 
sults in variations of 30-80 ms _1 (ij3] above), we instead use the 
relative tilt in each wavelength solution to correct the master wave- 
length frame (chosen as the wavelength solution at the start of ob- 
servations for each target star). Since the tilt is derived by com- 
bining all ThAr orders, we then apply this uniformly to each or- 
der in turn for the relevant ThAr observation that is made between 
each observation pair. This procedure is effectively equivalent to 
allowing only the low order terms to vary in the 2D wavelength so- 
lution. Each wavelength corrected frame is then used to make the 
subsequent deconvolved profiles (see next £14 . 1 1 ) for the pair of ob- 
servations that bracket it. The local shift for a given observation is 
made by measuring the shift of the telluric lines since this shift is 
not the same as that seen in the ThAr frames. Since this latter shift 
is assumed to apply to both the stellar and telluric lines, it is lim- 
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Figure 2. Stability of MIKE on the 21st (top) and 22nd (bottom) measured 
using cross-correlation of ThAr lines (crosses joined by solid line). Velocity 
in ms~ 1 is plotted against observed frame sequence number. Also plotted 
are the corresponding temperature variation of the dewar and the position 
(as altitude and azimuth) of the telescope. The temperature is measured in K 
and multiplied by 1000, while the azimuth can be read directly in degrees. 
The telescope altitude is multiplied by 4 for clarity. The hatched regions 
indicated periods over which the telescope was pointing at a fixed RA and 
Dec on the sky (i.e. tracking), and over which no slewing of the telescope 
took place. 

ited by the stability of the telluric lines. However, the correction, 
even for significant shifts only affects the zero point wavelength 
in our calculations. For example, even a shift of 1 pixel (equiva- 
lent to 2124 ms" 1 ) at 8000 A modifies the zero point wavelength 
by A A = XqAv/c = 0.057 A. The corresponding velocity error 
arising from a 0.057 A shift error (An = 8000 ± 0.057) is only 
~0.015ms _1 . 
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Figure 3. Comparison of Telluric shifts and ThAr shifts. The telluric lines 
positions were obtained from the second pair of observations of GJ 1061 
(our target with the highest S/N) on the first night. The shifts are plotted 
relative to the first observation of GJ 1061 as plus symbols (red) and crosses 
(green). Similarly the shift in the ThAr positions is plotted for the ThAr 
frame taken between the GJ 1061 observation pair relative to the first ThAr 
frame (taken after the first GJ 1061 frame) and are plotted as small points 
(blue). 

possesses a rich source of absorption lines. However, as Fig.[T]illus- 
trates, the I2 lines do not extend to the wavelengths we are using. In 
the near infrared, there has been some effort ma de towards identify- 
ing suitable cells (e.g. lMahadevan & Gel2009T) . while the only suc- 
cessful infrared gas cell s urvey has utilised a NH3 gas cell that pro- 
vides lines in the K band jBean et all20ich . The use of a laser comb 
to provide regularly spaced reference lines dSteinmetz et alj|2008h 
has now demonstrated 10 ms -1 on stellar targets (lYcas et al.l2012l) , 
although this technology is still being developed (at infrared wave- 
lengths) and currently remains an expensive choice. The only avail- 
able option for our data is to investigate the possibility of making 
use of the numerous telluric reference lines available in our spec- 
tra. Hence the intuitive requirement of stellar spectral regions that 
are void of atmospheric lines is replaced by a strategy that makes 
full use of the recorded spectral range in each exposure. We are 
therefore faced with a situation where the observed spectrum and 
reference lines are often superimposed, in much the same manner 
as when using a gas cell inserted into the instrument light path. Our 
adopted method is dictated by the nature of the spectra in the sense 
that we do not possess a reliable stellar template spectrum that is 
free of the reference spectrum (the telluric lines in this case). While 
accurate reference spectra can be derived from atmospheric mod- 
els, the stellar spectra, at least at high resolution, are not highly ac- 
curate, necessitating a semi-empirical approach. Rather than cross- 
correlating small regions of spectra and optimally co-adding, we 
employ a hybrid of the iodine t echnique of Marcy e t al. and the 
ThAr method first described by iBaranne et al. and now in- 

corporated into the HARPS pipeline. 



4 PRECISION RADIAL VELOCITIES OF M DWARF 
SPECTRA 

Since we are unable to achieve the desired sensitivity for preci- 
sion radial velocities using the internal ThAr emission line source, 
an alternative simultaneous fiducial is needed. To date, the only 
gas cell to be identified for precision work at optical wavel engths 
uses iodine (I2) as a reference fiducial dMarcv&Butlerlll992l) as it 



4.1 Least Squares Deconvolution 

We have applied least squares deconvolution (LSD) to our spectra 
to derive a single high S/N line profile for each observed spectrum. 
In order to disentangle the stellar signature from the telluric signa- 
ture, we apply the technique using models that represent the stellar 
lines and the telluric lines individually i.e. a stellar line list and 
telluric line list giving line positions and normalised line depths. 
Hence, two simultaneous high S/N ratio line profiles are derived, 
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one from the stellar lines and one from the telluric lines. The de- 
convolved profiles are derived in velocity space, enabling direct rel- 
ative shifts to be measured. 

The deconvolution procedure has bee n used in a number of 
applications but was first implemented by Don ati et"aT] i ll 9971) to 
boost the S/N in weak Stokes V profiles in noisy, high resolution 
timeseries spectra, thereby enabling indirect stellar magnetic im - 
ages to be derived. It was also implemented bv lBarnes et al. I dl998l) 
to similarly enable inversion of intensity profiles to produce high 
resolution surface brightness maps of relatively faint, rapidly rotat- 
ing Alpha Persei G dwarfs. More recently, the procedure has been 
used to search for the weak signature of planetary signals in high 
resolut ion optical ((Collier Cameron et al. 1999, 2002; Leig h et al.l 
l2003bl) and infrared dBarnes et al.ll2007aftbi 120081 12010|) spectra 
LSD is also effective for deriving accurate v sin i values for low 
S/N spectra that do not necessarily contain strong photospheric 
lines. The v sin i values of a nu mber of M dwarfs were derived in 
this way in ljenkins et al.l d2009h . We employ the same technique in 
Jenkins et al. (in prep.) to derive the v sin i values for our sample of 
M dwarfs in this paper. 

Although the method has been described before, here we give 
a detailed account of our implementation to emphasize its applica- 
tion to precision radial velocity measurements. If the deconvolved 
profile elements are denoted by Zk, then the predicted data, pj (i.e. 
the convolution of the line list and deconvolved profile), can be 
written as 



matrix, a, is effectively a weighting mask which only includes the 
regions of spectrum over the desired deconvolution range in wave- 
length/velocity space. 

Simply convolving the normalised, observed spectrum (ele- 
ments, rj) with the matrix, a, i.e. 



(5) 



would be equivalent to co-adding all the weighted lines in the object 
spectrum which occur in the synthetic line list, but would not prop- 
erly account for the effects of blended lines. Blending occurs in all 
spectra to some degree, being dependent on instrumental resolution 
and astrophysical properties, including stellar rotational line broad- 
ening. Blending occurs when there is more than one set of non-zero 
k elements per j element in the design matrix. The least squares de- 
convolution process means that the number of lines which can be 
used is independent of the degree of blending. 

The output least squares profile is binned at the average pixel 
resolution of the object spectrum and is determined by the size of 
the CCD pixels (approx. 2.1 -2.2 kms -1 in the case of our MIKE 
observations). Generally the velocity range over which the decon- 
volution takes place is chosen to include continuum on either side 
of the profile. The % 2 function is used to obtain the inverse variance 
weighted fit of the convolution of the line list with the deconvolved 
profile to the normalised and continuum-subtracted spectrum rj , 



Pi 



OljhZh, 



(1) 



The elements of the Nj x Nk matrix, a, are defined by the line 
depths, di, and the triangular interpolation function, A(x), as 



E\ 2 



(6) 



The equation is minimised by finding the solution to the set of equa- 
tions 



where 



Oijk 



x = \ v h - c 



X *~ Xi 1 1 /Av, 



(2) 



(3) 



and 

di = depth of line i 

Ai = wavelength of line i 

\j = wavelength of spectrum pixel j 

Vh = radial velocity of profile bin k 

Av = velocity increment per pixel in deconvolved profile 

The triangular function A(a;) is such that 

A(x) =l+x for -I<a:<0 
= 1 - x for <x <1 
= elsewhere 



(4) 



In this way, the weight of each spectrum element, j, is partitioned 
across 2 or more pixels, k. This interpolation in velocity space from 
each observed line to the deconvolved line thus takes account of 
the change in dispersion seen across the spectral orders. Hence a 
sparse matrix a is built up such that it contains minor diagonals, 
usually with pairs of partitioned elements in k for each element or 
pixel in j. In other words, the method given in equations |2]|3] and 
[4] is defined such that the interpolated pairs of elements in k, at 
a given j, add to give the depth of line i in question. The design 



dz k 



= 0, 



(7) 



where the unknown quantity is the least squares profile, Zk- This 
yields 



(8) 



The deco nvolution pro c ess c an be given in the same form as Equa- 
tion 4 in iDonati et al. I dl997h . If the inverse variances, 1/cr, are 
contained in the the vector, V, and the spectrum elements rj in 
the vector, R, the solution can be written in matrix form as 



z = {o^.V.a)- 1 a T .V.R 



(9) 



The right hand part of Equation |9](i.e. a T .V.R) is the cross- 
correlation of the observed spectrum with the line mask. This, as 
mentioned above, is the weighted mean of all the lines. The solution 
however gives flat continuum outside the least squares deconvolved 
profile, free of side lobes from blends, provided that the line list 
used is reasonably comprehensive, and the weights correct. The so- 
lution therefore effectively cleans the cross-correlation vector from 
the square symmetrical auto-correlation matrix a T .V.a. The so- 
lution to equation|9]requires the inverse of the auto-correlation ma- 
trix to be determined. This can be found by making use of a suitable 
fast inversion rou tine such as Cholesky Decomposition as given by 
|Pressetal.lfl986h . 
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Table 2. Radial velocities and cross-correlation uncertainties for each target. 



4.2 Deconvolution line lists 

For successful deconvolution, we require a good line list to 
represent both the synthetic telluric spectrum and the stellar 
spectrum. Since we find that spectra for mid-late M dwarfs 
dBrott & Hauschildtl [2005 ) do not adequately reproduce the line 
positions and strengths at the resolutions we are working at, we 
have adopted a semi-empirical approach for deriving line lists. We 
shifted and co-added spectra for one of our higher S/N targets (GJ 
1002). Once a high S/N ratio template is obtained (S/N = 320 for all 
co-added GJ 1002 frames), we simply search for absorption lines 
and find the line position by cubic interpolation of the line core. 
The line depths are measured at the same time for weighting in the 
deconvolution. We also included a S/N selection criterion to avoid 
selecting noise features as lines. Since we do not use lines that are 
weaker than 0.05 of the normalised depth, this should be avoided 
anyway since even lines with a depth of 0.05 are 16 sigma above the 
noise level. Finally, the lines in the stellar line list that were either 
close to known telluric lines (see below), or actual telluric lines, 
were rejected to remove the possibility of cross-contamination. The 
line list is not ideal since it is derived at the same resolution as the 
observation, whereas a higher resolution spectrum would be desir- 
able to disentangle blends. Nevertheless, any incorrect wavelength 
positions from blended lines will be treated the same way in every 
deconvolution of a given target star. 

To generate the synthetic telluric spectra, we use the Line 
By L i ne Ra diative Transfer Model (LBLRTM) code dClough et~al] 
1 1991 12005b . which is distributed by Atmospheric and Environ- 
mental Research (AERfj. The code uses the most up to date 



versions of the HITR AN 2008 molecular spectroscopic database 
dRothman et~ai] ?2009). We generate a reference spectrum at high 
resolution, using a meteorological sounding provided by Air Re- 
sources Laboratory (ARLjfl The soundings are available at a time 
resolution of 3 hours for the longitude and latitude of Las Cam- 
panas Observatory. The sounding contains an atmospheric obser- 
vation containing wind speed, wind direction, temperature, dew 
point and pressure at a range of atmospheric heights and is used 
by LBLRTM to generate synthetic spectra for the current location 
and conditions. Our telluric line list is generated from a reference 
synthetic spectrum with a sounding taken on the first night of obser- 
vations at Las Campanas observatory for an altitude of 60 degrees. 
Once a high resolution normalised spectrum is generated, we build 
a line list in the same manner as described above for the stellar 
lines, although contamination and S/N effects in this instance are 
not relevant. 



4.3 Derivation of radial velocities 

Radial velocity measurements are made by subtraction of the mea- 
sured velocity centre of the deconvolved stellar profile with respect 
to the deconvolved telluric profile for each s pectru m. We use a ver- 
sion of the HCROSS a lgorithm of Heavens ( 1993) which is in turn 
a modification of the iTonrv & Davisl ( 19791) cross-correlation al- 
gorithm. HCROSS applies a more accurate, robust algorithm that 
employs the theory of peaks in Gaussian noise to determine uncer- 
tainties in the cross-correlation peak. Since the routine, which is 
part of the Starlink package , FIGARO (now distributed by the Joint 



2 http://rtweb.aer.com 3 http://ready.arl.noaa.gov/READYamnet.php 

© 2010 RAS, MNRAS OOO.fTlHll 



8 J.R. Barnes et al. 



GJ 1002 Tellurics 



E 




I I 0.2 I I 

-30 -20 -10 10 20 30 -30 -20 -10 10 20 30 

Velocity [kms -1 ] Velocity [kms -1 ] 



Figure 4. Right panel: Deconvolved profiles for GJ 1002 plotted in order 
of observation from top to bottom. The left panel shows the simultaneously 
deconvolved telluric profiles. Error bars are also plotted but too small to be 
visible. 

Astronomy Central), was originally intended for measurement of 
galaxy redshifts in low S/N spectra, we have made a minor modifi- 
cation to output both shift, and shift uncertainty, in pixels. 

A master stellar line is derived by cross-correlating and align- 
ing all deconvolved stellar lines. The master line is used for as a 
cross-correlation template since it more closely represents each in- 
dividual line as opposed to using a simple Gaussian, Lorentzian 
or Voigt profile. A similar procedure is carried out for the telluric 
lines before the subtraction of pairs of stellar and telluric measure- 
ments can be made. Finally, we apply the bary centric corrections 
at the mid-time of each observation to correct our velocities to the 
solar-system centre. The final velocities are arbitrary and depend on 
the stellar line lists used (i.e. the radial velocity of GJ 1002). With 
our small sample spanning two days, we have opted to subtract the 
mean velocity from all data points of each target. With higher reso- 
lution template spectra, we will determine the absolute barycentric 
radial velocity of each observation in future work. 



5 RESULTS 

5.1 Target star radial velocities 

Our small sample of 7 bright mid-late M dwarf stars enabled one 
to two sets of observations to be made on each night. An exam- 
ple of the deconvolved profiles, plotted for GJ 1002, is shown in 
Fig- 111 To monitor stability and reproducibility of results, we ob- 
served each target twice in a single pointing, with an intermediate 
ThAr frame, as discussed above. After extraction, we obtained S/N 
ratios of ~25 - 160, which after deconvolution yielded single lines 
with S/N ratios of 875 - 7400. The vital statistics of each star, in- 
cluding the total number of observations are listed in Table [T] The 
radial velocities derived from our targets are listed in Tableland 
plotted in Fig. [5] Statistics, pertaining to precision of the measure- 
ments are given in Table[3] We are confident that the error estimates 
from HCROSS are a good representation of the cross-correlation 
uncertainties. This is illustrated by columns 4 and 5 which give 

4 http ://starlink.j ach. hawaii . edu/ starlink 



the HCROSS error and the mean difference between pairs of obser- 
vations for each object. In other words, over the short timescales 
on which pairs of observations are made (with no slewing of the 
telescope), the scatter in the points agrees well with the cross- 
correlation error estimate. While this also holds for GJ 3128 when 
the last observation pair (with a difference of 42 ms _1 ) is excluded, 
the much higher v sin i and low S/N ratio of LP 944-20 may be the 
cause of the breakdown of this one-to-one relationship. There are 
also only two pairs of observations of LP 944-20 since the first was 
a single observation. 

In Table [3] the final column lists the discrepancy between the 
HCROSS error and the scatter. The value is simply that which added 
in quadrature to the HCROSS error leads to the observed scatter in 
the observations. As such it represents any further unaccounted for 
errors in the measurement, which may be astrophysical, instrumen- 
tal, a real planetary signal, or a combination of these factors. In the 
next section we discuss such possible causes further. 

The flattest RV curve is exhibited by GJ 1061, with an overall 
r.m.s. scatter of 15.7 ms _1 , comparable with the cross-correlation 
uncertainty of 9.37 ms -1 derived via HCROSS. While the GJ 1286 
RV curve is relatively flat, the overall scatter is more than twice 
the estimated error. In the case of GJ 1002, this discrepancy is 
closer to a factor of 3. Although LHS 132 appears relatively flat, 
we again emphasise that there are only two pairs of observations 
of this object, both taken at almost the same time on each night. 
No discernible variation is found in LP944-20 which shows the 
largest scatter and errors, owing to its late spectral type and high 
usini. To the contrary, GJ 3128 shows the same rising trend on 
each night of observations with a scatter that is ~10 times the er- 
ror estimate. Perhaps the most interesting RV curve is exhibited by 
SO J025300+165258, where a rising trend is seen on night one fol- 
lowed by a falling trend on the second night. We discuss GJ 3128 
and SO J025300+165258 in more detail in section ( jj6T) . 

5.2 Photon noise sensitivities 

To validate the technique we have carried out Monte-Carlo sim- 
ulations to determine the radial velocity precision we expect to 
achieve as a function of S/N ratio. Synthetic spectra at R = 35,000 
were generated from the semi-empirical line lists that we derived 
in i]4.2l The spectra contained both telluric lines and spectral lines. 
An observed continuum (of GJ1002) was used to mimic the counts 
across the whole spectrum, thereby accounting for the spectral type 
and instrumental response (e.g. quantum efficiency and blaze func- 
tion) of MIKEat R ~35,000. Noise was then added to the spectrum, 
and subjected to out pipeline procedure for deriving radial veloc- 
ities. This procedure was carried out 100 times at each S/N level, 
resulting in 100 radial velocities from which the scatter, indicating 
the measurement precision was derived. The results are shown for 
mid-M dwarfs in Fig.|6]for non-rotating stars. For the S/N range of 
25 - 160 seen in our sample of stars, the photon limited precision is 
1.7-11 ms _1 respectively. The right panel of Fig. [6] demonstrates 
that with 10 ms~ 1 precision (for example by tellurics), there is little 
benefit in obtaining spectra with S/N ratios > 100. 

The lower curves plotted in Fig. [6] (line joined by open cir- 
cles) show the precision attainable from tellurics alone, indicat- 
ing that they possess the larger contribution for low v sin i (i.e 
ii sin i ^ 5 kms -1 ). For zero rotation (vsini = 0kms _1 ), the tel- 
luric cross-correlation uncertainties in the simulation are ~ 1.5 
times the stellar uncertainties. For example, with S/N = 120 and 
i;sini = 0kms~ , the telluric and stellar contributions are 1.5 
ms _1 and 2.2 ms _1 respectively. The combined precision of 
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Figure 5. 2010 November 21 & 22 heliocentric corrected radial velocity measurements for the 7 M dwarfs, GJ 1002, GJ 1286 and GJ 1061, GJ 3128, SO 
J025300+165258, LHS 132 and LP944-20. The measurements were made by deconvolving the spectral lines in each frame into a single high S/N profile. 
The same procedure is applied to each telluric frame. The horizontal dashed line simply represents the mean velocity (zero). The pre-deconvolution S/N ratio, 
mean cross-correlation error, cr, and r.m.s. scatter in all points are given for each star. For GJ 3128 and SO J025300+165258, which both exhibit significant 
radial velocity variations, we fit a sinusoidal velocity curve (see §6.1 for discussion). 
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2.6 ms -1 is similar to that achieved using an I2 cell and VLT-UT2 + 
UVES for the Barnard 's star (M4V ) (Kiirster et al. 2003) and Prox- 
ima Centauri (M5V) dEndl & Kurstedl2008l) . We note however that 
both objects possess earlier spectral types and additionally, were 
observed using an image sheer at a much higher resolution of R ~ 
100,000. A direct comparison of precision and efficiency (see also 
S]2.2b requires further simulations at other wavelengths to be made, 
although increased resolution would clearly improve our simulated 
precision further. We reiterate that telluric stability will neverthe- 
less ultimately limit precision to a few ms - . 



6 DISCUSSION 

While precision of a few 10s of ms - has been demonstrated by 
our results presented in Fig. [5] it is clear that the observations do 
not reach the Poisson noise limit. Our simulations (Fig. [6) indicate 
a photon noise limit of order 2.4 ms~ for our brightest target, GJ 
1061 (i.e. at S/N = 163). For GJ 1002 (S/N =113), the photon noise 
limit is 3.3 ms" 1 and similarly for GJ 1286 and GJ 3128 (S/N = 
84), the limit is 4.4 ms" 1 . For a modest S/N = 25, as for LHS 132, 
we expect a limiting precision of 14.5 ms -1 . With the exception 
of this latter case, for which we anyhow only have two observa- 
tions, the photon-noise limit is clearly not sufficient to explain the 
discrepancy between the r.m.s. scatter in the points and the error 
bars listed in Table |3] LP944-20 is a rapid rotator with v sin i ~ 35 
kms -1 and our estimated photon noise limit (not shown in Table 
at S/N = 44 is 99 ms~ and therefore in reasonable agreement 
with both the error and scatter, which is smaller than the difference 
between pairs of observations. More data is needed on this kind of 
object whic h has a high v sin i(Jenkins et al., In prep.), even among 
M dwarfs jjenkins et alj2009h . 

The question of telluric stability has been tackled in a num- 
ber of publications. Measurements of Arcturus and Procyon w ere 
made with ~50 ms" 1 precis ion by iGriffin & Griffinl dl973l) . A 
number of authors, including Snellen! ( 120041) have demonstrated 
that sub-10 ms -1 is achiev able by observin g H2O lines at wave- 
lengths longer than 0.6 urn. iGrav & Brownl d2006l) however found 
25 ms -1 precision from tel luric observations of r Ceti. More re- 
cently Fig ueira et alj d2010i) have investigated the use of telluric 
features over timescales of a few days to several years by measuring 
the stability of some of the stronger O2 lines found in HARPS obser- 
vations. Their observations indicate that for bright targets such as 
Tau Ceti, precision of ~ 5 ms -1 is achievable on timescales of 1 - 8 
d when observing at airmasses restricted to < 1.5 (altitude >42°). 
Inclusion of data at all air masses (<2.2, altitude > 27°) enables 
precision of ~ 10 ms -1 to be achieved. Over longer timescales of 
1-6 years, the same precision of ~ 10 ms -1 is found. In other 
words, atmospheric phenomena are found to induce radial velocity 
variations in telluric lines at the 1 - 10 ms -1 level, but can largely 
be corrected for. Moreover, and importantly, simple fitting of at- 
mospheric effects, including asymmetries found in molecular lines 
that vary as a function of altitude, and wind speed, enable correc- 
tions down to ~ 2 ms -1 to be made. Figueira et al. find that this 
is twice the photon noise limit of their observations. If we include 
a 10 ms -1 uncertainty into our error budget, our overall contribu- 
tion from tellurics and photon noise ranges from 10.2 ms -1 (GJ 
1002) to 17.6 ms -1 (LHS 132). This simple argument is almost 
sufficient to declare the radial velocity curve of GJ 1061 to be flat 
within our uncertainties since the estimated discrepancy listed in 
Table[3]is 12.6 ms" 1 . Again, the same argument may be made for 
LHS 132, but clearly more observations are needed. The remaining 
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Figure 6. Left: Simulation showing the Poisson noise limited precision of 
our radial velocity method applied to MIKE for 0.62-0.90 urn with R = 
35,000. The simulation includes both the stellar and telluric line contri- 
butions. The curves are plotted for stellar rotation of usinz = 0, 5, 10 & 
20 kms" 1 . Right: The same curves showing the simulation with a 10 ms 
telluric uncertainty added in quadrature. In both panels the lowest curve 
(open circles) indicates the limiting precision which is set by the telluric 
lines. 



M5.5V stars require modest (GI 1286) to significant (GI 1002) ad- 
ditional uncertainties to explain their discrepancies, while GI 3128 
and SO 1025300+165258 both show more significant discrepan- 
cies. It is additionally worth noting that all observations of GJ 3128 
were observed in the airmass range 1.38 -1.43 (44° - 46°), while SO 
1025300+ 165258 was observed with an airmass range of 1.49-1.71 
(36° -42°), so airmass effects are unlikely to lead to radial velocity 
variations of order ~ 100 ms - (see i j6.1l l. 

To investigate further the expected errors that may arise from 
molecular asymmetries in telluric lines, we used LBLTRM to simu- 
late telluric spectra at an airmass range of 1.0-2.2. At R = 35,000, 
we find only ~ 1 ms -1 variation for observations made at Las Cam- 
panas. The result additionally gives a measure of any uncertainty 
from using a mismatched line list for deconvolution. All our re- 
sults have used a master telluric line list for an altitude 60° . Hence 
we do not expect airmass corrections to be important until we are 
confident that we have reached precision of a few ms -1 . A more 
important consideration potentially arises from wind speed and di- 
rection. Since most of the lines in our line list are water transitions, 
we expect the lines to form at lower altitudes. The soundings used 
(see &l4.2t by LBLTRM give wind direction and speed for altitudes 
in the range ~2 - 27 kms -1 . Since the water lines form largely at 
lower altitude, we note that the wind speed over the two nights 
of observations is < 10 ms -1 on both nights. Additionally, the 
wind direction changes by only ~20°, giving a maximum differ- 
ence between the two nights of ~ 0.6 ms" 1 . With a wind vector of 
360° (Northerly), at altitudes <5000 m, the maximum wind speed 
seen by a star at altitude 40° in a 10 ms -1 wind is 10cos(40°) = 
7.7 ms -1 during our two night observing run. Of course in real- 
ity, the differential is somewhat less since high declination stars 
only traverse a small range of azimuth angles. For a northerly wind 
at Southern latitudes, intermediate latitude stars traverse a greater 
range of azimuth angles, but will never see as much as the predicted 
full 7.7 ms -1 variation. We are therefore confident that in the ab- 
sence of high precision measurements, telluric velocity corrections 
can not give the radial velocity variations seen in some of our stars. 
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Object S/N Error Mean AOP r.m.s. scatter Discrepancy Photon limited 

[ms _1 ] [ms _1 ] [ms _1 ] [ms _1 ] precision [ms _1 ] 



GJ 1002 


113 


12.8 


12.7 


32.4 


29.7 


3.3 


GJ 1061 


163 


9.37 


11.8 


15.7 


12.6 


2.5 


GJ 1286 


84 


8.63 


8.88 


22.1 


20.3 


4.4 


GJ3128 


84 


8.80 


20.0(12.6) 


87.7 


87.2 


4.4 


SO J025300+165258 


90 


9.91 


7.0 


63.3 


62.5 


4.2 


LHS 132 


25 


29.7 


27.2 


33.7 


16.0 


14.5 


LP944-20 


44 


69.1 


208 


121 


99.3 


7.6 



Table 3. Radial velocity errors for each target. The third column is the error derived from the cross-correlation analysis. Column 4 (AOP) is the mean 
uncertainty of the pairs of radial velocity points. Column 5 is the r.m.s. scatter of all observations for the object indicated. The value in brackets for GJ 3128 
excludes the final pair which result in the higher value of 20.0 when all pairs are considered (see also Fig. 4). Column 6 lists the discrepancy between the 
scatter and the errors (columns 3 and 4), indicating the magnitude of other radial velocity components (i.e. planetary signal or noise) and column 7 gives the 
photon limited precision with MIKE. 



As MIKE does not have an image derotator it is not able to keep 
the slit at the parallactic angle (it is set such that the slit is vertical 
at 30° from the Zenith (altitude = 60°). It is unclear whether this 
should have an effect on the measured radial velocities. One might 
reasonably expect that to first order, any slit effects would cancel 
owing to the reference lines and stellar lines following the same op- 
tical path in the instrument. There is no clear systematic correction 
that can be applied to all the objects that appears to correct, in an 
empirical manner, for the slit tilt. 



6.1 Planetary candidates? 

Although GJ 3128 appears to show significant radial velocity vari- 
ations that are at 10rr of the cross-correlation errors, the same as- 
cending RVs of both observation pairs on both nights raises the 
possibility of an as yet unidentified one day aliasing effect. The 
radial velocity points can be fit with a simple sinusoidal curve 
with a period of P = 0.511 d, an amplitude of K* = 209.4 ms _1 
and a residual r.m.s. of 1.8 ms _1 . A circular (eccentricity, e = 0) 
planetary orbit of this period and amplitude could be induced by 
a m p sin i = 3.29 Mn cp planet in a very close 0.0058 AU orbit 
(i.e. 8.31 R*). We note also that under the assumption of a max- 
imum v sin i = 5 kms -1 , the minimum period of the GJ 3128, wit h 
R = 0.15 R dReid & Hawievll2005l ; iKaltenegger & Traubll2009h . 
would be P ~ 1 .52 d, quite distinct from the RV signal. Importantly, 
Fig. [7] shows that GJ 3128 exhibits flaring activity. Our second ob- 
servation shows clear evidence of such an event which has pushed 
Ha int o emission. Thi s interesting phenomenon confirms the find- 
ings of IReinersl d2009h . already discussed in §1, and demonstrates 
that we are not sensitive to flaring events at the level of our pre- 
cision since the first pair of points do not show significant radial 
velocity differences. Moreover, since no further strong flaring is 
seen in the remaining spectra, we believe that such events should 
not be responsible for the RV variation seen in GJ 3128. We also 
emphasise that we exclude a region of ± 7 A surrounding Ha from 
our deconvolution. 

Our most promising candidate for a stellar radial velocity 
signature is exhibited by SO J025300.5+165258, which shows 
a rising trend on the first night and a falling trend on the sec- 
ond night. The r.m.s. scatter is at a level of 6.4<r of the cross- 
correlation errors. A sinusoidal fit to the data indicates a period 
of P = 2.06 d with an amplitude of K, = 182.7 kms" 1 (residual 
r.m.s. = 0.46 ms -1 ). A circular planetary orbit would lead to a 
mass of m p sin i = 4.88 Mn cp with a = 0.014 AU, falling inside the 
classical habitable zone. We emphasise that we make no claims 
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Figure 7. The normalised spectral region around Ha for all GJ 3128 and 
SO J025300.5+165258 spectra. For GJ 3128, the second observation from 
the first night shows clear evidence of a flaring event which has pushed Ha 
into emission. SO J025300. 5+165258 shows that Ha is variable from one 
spectrum to the next, over all 9 spectra. No other spectra show this degree 
of variability. 



for a planet and clearly further observations are required. In fact, 
perhaps Fig. |7] shows that in fact SO J025300.5+ 165258 is the 
most chromospherically active sta r in ou r sample, but on the ba- 
sis of our GJ 3128, the IReinersl j2009h results and removal of 
these regions from the RV analysis, it seems unlikely that these 
chromospheric variations are responsible for the RV variations of 
SO J025300.5+165258. However, while pairs or triplets of observa- 
tions are consistent at the cross-correlation level, the Ha line also 
does not show large variations on these timescales, but does in fact 
vary from one pair/triplet to the next. Additionally, for complete- 
ness, we note that no other stars exhibit significant variations in the 
Ha line region, except for the final pair of observations of GJ 1286, 
which show a slight increase Ha emission. 
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6.2 Line bisectors 

Although we believe that chromospherically active lines them- 
selves do not bias our radial velocities, the active nature of the star 
clearly warrants further analysis of activity indicators. Line bisec- 
tors are a commonly employed tool that enable line asymmetries 
due to photospheric temper ature inhomogeneities such as starspots 
and/or plage to be assessed ( Saar & Donahue 1997). Following the 
same procedure as lMartfnez Fiorenzano et al.l (120051) . we calculate 
line bisectors to assess starspot jitter. Fig. [8]shows our results for 
SO J025300.5+165258 for both the deconvolved stellar lines and 
telluric lines. The line bisectors are shown for all 9 profiles. The 
bisector span is calculated by determining the mean bisector for 
two regions of the line profiles, near the top and the bottom of the 
profile. Thus, any spot or spot group on the star that rotates into 
and out of view will cause variation in the span, with the subse- 
quent line asymmetries yielding false RV signals (as defi ned, the 
bisec tor span, BIS, shows an anti-correlation with RV (Des ort et al.l 
2007)). The telluric line profiles are quite asymmetric in the line 
wings, but are relatively consistent over all observations. The stel- 
lar lines show more variation in the wings (i.e. at normalised depths 
>0.6), but even points 5-7 (the triplet of observations on night 2 
which show bisector morphology changes in short timescales), do 
not show large variation in their RVs. 

With so few observations, we find that it is in fact not possi- 
ble to find a bisector-span vs RV correlation owing to the scatter 
in the bisector values. We combine bisector span values from both 
the telluric and stellar lines since this should also account for any 
changes in the telluric lines. With a simple linear relationship, it is 
possible to flatten the two pairs of points on night 1, but this has lit- 
tle effect on the points from night 2. It is clear that a larger number 
of data points are needed at more phases before any bisector span 
vs radial velocity relationship can provide useful information on 
the nature of radial velocity variations. However, with a total varia- 
tion in the bisector span of ~150 ms _1 for SO J025300.5+ 165258 
(92 ms -1 when the large span calculate for profile 6 is removed), 
we can't completely rule out that a trend would emerge with more 
observations. We also find similar bisector variations in our other 
stars, including our brightest target, GJ1061, which is nevertheless 
relatively flat within the estimated cross-correlation uncertainties. 
It is worth noting that the ~1 per cent FWHM variations that we 
see in our profiles also do no t yield a clear (anti-)correlation, as re- 
ported bv lBoisse et al] < |201 lh . As with SO J025300.5+165258, we 
are unable to determine a correlation for GJ 3128, but do find that 
the second observation (where we reported a flare, as shown if Fig. 
[7), is consistent with the first observation, confirming that it had 
little effect on the photospheric lines. The above findings may well 
be a consequence of lower contrast starspots in low temperature at- 
mospheres that contribute negligible stellar jitter, a hypothesis that 
may yield clearer answers with further data. 



6.3 Detection analysis and limits 

In the case of no planetary RV variations, we are sensitive (at la 
scatter) down to plan et masses of 6 in the case of GJ 1061. 
Previous simulations dBarnes et alj|201 lh in the near infrared Y- 
band have shown that only of order 20 epochs are required to de- 
tect 5 M9 planets orbiting 0. 1 Mq stars rotating with v sin i= 5 
kms -1 a nd exhibiting solar- like spot activity. The simulations pre- 
sented in lBarnes et al.l d201 lh used standard periodogram analyses 
to search for planetary signals in fake radial velocity curves injected 
with astrophysical starspot jitter. However an analysis of several 
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Figure 8. Line bisector plots for all 9 observations of SO 
J025300.5+165258 for the deconvolved stellar (left) lines and telluric 
lines (right). We are unable to determine a clear bisector span vs RV 
correlation, but note that we observe 143 ms variation in the stellar 
bisector span. 

marginal detections and non-detections using Bayesian techniques 
that are efficient at searching for low amplitude signals in nois y data 
dKotiranta & Tuomj|20ld : lTuomil201 ll i lTuomi & Jonesl201 ll) indi- 
cate that the number of epochs required is conservative. Clustering 
of randomly located starspots in our 3D stellar models was found 
to lead to non-symmetric noise/jitter distributions using these rou- 
tines. In some instances, using the correct noise-model resulted in 
the recovery of orbital parameters, where no planet had even been 
detected at the 1 per cent false alarm probability level of our more 
straightforward periodogram analysis. This observation may well 
prove vital for stellar systems where the rotation period of the star 
may potentially be close to that of the planet, a distinct possibility 
for HZ planets orbiting mid-late M stars. 



7 SUMMARY & PROSPECTS 

By making use of the red optical spectral regions and an efficient 
CCD detector we have demonstrated that sub-20 ms~ precision 
radial velocity measurements of M dwarfs are possible when utilis- 
ing a telluric reference fiducial and a wavelength extent (~3000 A) 
equal to more traditional optical surveys. Two stars show signifi- 
cant radial velocity variations that are consistent with the RV am- 
plitudes that could be induced by close orbiting planets of several 
Neptune masses. At this stage, with so few observation epochs, 
we do not argue that these are planetary in origin, although we 
note that Neptune mass planets are expecte d (Howard et al .1 120 101 ; 
IWittenmver et al.l 1201 ll ; iBonfils etallboill ; iBorucki et aimoill) 
with a frequency that is consistent with detection rates of 1/7 - 2/7. 
While we find that the chromospheric activity, seen as variation in 
Ha emission, is unlikely to be the cause of the variations, based 
on local consistency of observation pairs (e.g. where one observa- 
tion of a pair contains clear evidence for a strong flare) and other 
findings dReinersl2009l) . further observations are required to enable 
a fuller investigation aimed at the role of photospheric line shape 
variability. 

Using simulations based on our MIKE observations, we have 
shown that sub- 10 ms -1 precision can potentially be achieved for 
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spectra with S/N ~50. This improves to precisions of order 3 
ms _1 for S/N ~ 100 and 2 ms _1 for S/N ~ 200. Since rotation 
is known to increase on average for mid-late M dwarfs, more mod- 
est precisions of 5 - 12 ms~ are more likely once rotation speeds 
reach typical 10-20 ms _1 for M6V and M9V stars respectively 
dJenkins et ai]|2009l) . 

Our observations suggest that we have achieved reference 
fiducial limited observations of at least one target, GJ 1061, where 
the discrepancy between our cross-correlation errors, which are 
accurate on short timescales, and the overall scatter, can be ex- 
plained by atmospheric variations from tellurics. Since these should 
in fact be sub-10 ms _1 during our observations, there are likely ad- 
ditional uncertainties which most probably arise from MIKE. We 
have demonstrated that the instrument does not show a slow veloc- 
ity drift with time, but rather that apparent drifts of a few hundred 
ms _1 can be observed on relatively short timescales of a few ob- 
servations. We should add that MIKE was not built with precision 
radial velocities in mind and is not pressure or temperature sta- 
bilised, although the Magellan Planet Search program which oper- 
ates in th e optical and utilise s an iodine cell is a chieving 5 ms _1 
precision (Min niti et al.ll2009l ; lArriagada et alj2010h . It is not clear 
whether a mechanical effect, possibly introduced by the calibration 
mirror while taking arcs, or whether telescope slewing and possible 
disturbance of the Nasmyth mounted instrument are responsible for 
these shifts. Additionally, we are unable to determine whether the 
shifts are pseudo-random and occur only during movement of the 
telescope and/or instrument components, or whether the influence 
of slow temperature and pressure changes throughout the night are 
important. If gravity settling of the instrument is significant (esp. 
after a CCD dewar refill), we would expect possible drifts of up 
to several 100 ms _1 during a single observation. Fig. [2] appears 
to indicate that the shifts are not entirely random as they often con- 
tinue in one direction over several observations. As discussed in the 
previous section for the variable slit angle, owing to the simultane- 
ous measurement of reference fiducial and stellar spectrum, which 
traverse the same optical light path, such effects may well be ex- 
pected to cancel to first order. Nevertheless, they may go some way 
to explaining the discrepancy between residuals and our error bar 
estimates. 

In light of such difficulties, MIKE, operating at a modest res- 
olution of ~35,000 (especially when we consider that the tem- 
perature and pressure stabilised HARPS operates at R = 110,000) 
has demonstrated that precision radial velocities are well within its 
grasp. This work would nevertheless clearly benefit from a more 
stable instrument operating at higher resolution. We expect our pre- 
cision to scale approximately linearly with resolution, so the above 
estimates of photon limited noise, would potentially enable 1 ms~ 
precision to be achieved on the brightest slowly rotating M dwarfs 
at a resolution of R ~50,000. The additional S/N loss at this resolu- 
tion could easily be offset by telescopes with a larger aperture than 
the 6.5m Magellan Telescopes could offer. While we have not been 
able to determine the effects of a slit that does not maintain paral- 
lactic angle on the sky, we expect that doing so could only add to 
the precision (even if reddening effects at low airmass are reduced 
in the re-optical). Instruments such as UVES have already demon- 
strated precision of a 2 - 2 .5 ms~ over 7 years using an iodine cell 
dZechmeister et al]|2009h and additionally offer greater red sensi- 
tivity. This opens up the potential of using the 0.9- 1.0 \im wave- 
length region. The richness of atmospheric lines over most of this 
range, combined with the relatively free z-band region just short of 
1 .0 \im, and where M dwarf fluxes are even greater, would clearly 
be of great benefit to this kind of study. 



We see no rea son why the 2 ms _1 precision reported by 
Figue iraet al.l J2010h can not be achieved on existing 8 m class tele- 
scopes with existing instruments. By spectral type M6V, habitable 
zone rocky planets of 1 Me would be expected to induce radial 
velocity variations of similar magnitude to this precision while the 
most massive rocky planets with 10 would be capable of de- 
tecting signals at levels of order 10ct in ^ 20 epochs of observations 
over short timescales of days. 
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